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ABSTRACT 

We present the first hydrodynamic, multi-dimensional simulations of He-shell flash convection. 
Specifically, we investigate the properties of shell convection at a time immediately before the 
He-luminosity peak during the f 5 th thermal pulse of a stellar evolution track with initially two 
solar masses and metallicity Z = 0.01. This choice is a representative example of a low-mass 
asymptotic giant branch thermal pulse. We construct the initial vertical stratification with a set of 
polytropes to resemble the stellar evolution structure. Convection is driven by a constant volume 
heating in a thin layer at the bottom of the unstable layer. We calculate a grid of 2D simulations 
with different resolutions and heating rates. Our set of simulations includes one low-resolution 3D 
run. The computational domain includes 11.4 pressure scale heights. He-shell flash convection 
is dominated by large convective cells that are centered in the lower half of the convection 
zone. Convective rolls have an almost circular appearance because focusing mechanisms exist 
in the form of the density stratification for downdrafts and the heating of localized eddies that 
generate upflows. Nevertheless, downdrafts appear to be somewhat more focused. The He-shell 
flash convection generates a rich spectrum of gravity waves in both stable layers above and 
beneath the convective shell. The magnitude of the convective velocities from our ID mixing- 
length theory model and the rms-averaged vertical velocities from the hydrodynamic model are 
consistent within a factor of a few. However, the velocity profile in the hydrodynamic simulation 
is more asymmetric, and decays exponentially inside the convection zone. An analysis of the 
oscillation modes shows that both g-modes and convective motions cross the formal convective 
boundaries, which leads to mixing across the boundaries. Our resolution study shows consistent 
flow structures among the higher resolution runs, and we see indications for convergence of the 
vertical velocity profile inside the convection zone for the highest resolution simulations. Many 
of the convective properties, in particular the exponential decay of the velocities, depend only 
weakly on the heating rate. However, the amplitudes of the gravity waves increase with both the 
heating rate and the resolution. 

Subject headings: 



1. Introduction 

We present hydrodynamic simulations of the in- 
terior shell convection zone driven by the He-shell 
flash in thermal pulse Asymptotic Giant Branch 
(AGB) stars. 
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1.1. AGB evolution 

AGB stars are the final evolution stage of low- 
and intermediate mass stars before formation of 
a white dwarf (Iben & Renzini 1983). Nuclear 
production in AGB stars contribute to the chem- 
ical evolution of galaxies (for example Travaglio 
ct al. 2004; Renda et al. 2004). The chemical 
yields are support the interpretation stellar abun- 
dance from dwarf-spheroidal galaxy satellites of 
the Milky Way relative to stellar abundances of 
galactic halo stars (Venn et al. 2004; Geisler et al. 
2005). This may provide important clues about 
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the cosmological origin of our galaxy. Another 
role of AGB stars have recently been discussed as 
a potential source of the abundance anomalies ob- 
served in globular cluster star members (Ventura 
et al. 2002; Denissenkov & Herwig 2003), a ques- 
tion that is still debated. 

AGB stars produce substantial amounts of Li, 
C, N, 22 Ne, 23 Na, and the neutron-heavy Mg iso- 
topes. These stars are the progenitors of planetary 
nebulae nuclei and white dwarfs, including those 
thought to orbit C-rich extremely metal-poor 
stars with s-process signature (CEMP-s, Beers 
& Christlicb 2005). The favoured interpretation 
of their strongly non-solar abundance pattern is 
the pollution of the EMP star with material trans- 
fered from the white dwarf progenitor AGB star. 
This interpretation is largely based on the fact 
that the abundance patterns of CEMP-s stars in 
general resemble the C- and s-process overabun- 
dance predicted for EMP AGB stars. However, 
very significant discrepancies between observed 
and predicted abundances exist which could indi- 
cate some serious problem with current ID stellar 
evolution and nucleosynthesis simulations. Most 
likely those problems are related to the ID ap- 
proximate treatment of mixing that originate from 
multi-dimensional, convective fluid motions. 

Especially at extremely low mctallicity, AGB 
stars can produce many species in a primary mode, 
i.e. only with the H and He available from the Big 
Bang (Siess et al. 2002; Herwig 2004, for example). 
AGB stars are also the nuclear production site of 
the main component of the s process (Busso et al. 
1999), which produces roughly half of all trans- 
iron elements (Arlandini et al. 1999). The in- 
terpretation of high-precision laboratory measure- 
ments of isotopic ratios of heavy elements in pre- 
solar meteoritic SiC grains depends sensitively on 
the thermodynamic and the closely related mixing 
conditions at the bottom of the He-shell flash con- 
vection zone (Zinner 1998; Lugaro et al. 2003a, b). 
In this particular layer the 22 Ne(a, n) 25 Mg reac- 
tion releases neutrons very rapidly, and activates 
many s-process branchings. Experimental mea- 
surements of nuclear properties of the unstable 
branch-point nuclei are now underway (Reifarth 
et al. 2005) and in conjunction with grain mea- 
surements will provide new constraints for mixing 
at convective boundaries in the future. 

Currently, only ID stellar evolution models ex- 



ist to describe the processes of the AGB stellar 
interior quantitatively However, while spherical 
symmetry is globally a very good approximation, 
a detailed description of mixing processes induced 
by the hydrodynamic fluid-flows requires multi- 
dimensional simulations. Only over long times, 
equivalent to many convective turn-overs, can they 
be described accurately by the average quantities 
provided by the commonly used local theories, like 
the mixing- length theory (Bohm-Vitense 1958). 
However, it is well known that these theories do 
not describe the quantitative properties of mix- 
ing at convective boundaries. For these reasons, 
we start our investigations of the hydrodynamic 
properties of stellar interior convection with the 
He-shell flash convection in thermal pulse AGB 
stars. 

1.2. Hydrodynamic simulations of stellar 
convection and application to ID stel- 
lar evolution calculations 

Most multi-dimensional simulations of stellar 
convection address the outer convection zones. 
For example the solar convection zone includ- 
ing the surface granulation (Nordlund 1982; Stein 
& Nordlund 1998), convection below the surface 
(Stein & Nordlund 1989) or the role of mag- 
netic fields and rotation in the solar convection 
zone (Brun 2004) Other simulations include work 
on deep envelope convection in giants (Porter & 
Woodward 1994, 2000; Freytag et al. 2002; Robin- 
son et al. 2004), shallow surface convection, for ex- 
ample in A stars and white dwarfs (Freytag et al. 
1996) and the general properties of slab convec- 
tion (Hurlburt et al. 1986). Hydrodynamic sim- 
ulations of stellar interiors include an investiga- 
tion of semi-convection in massive stars by Merry- 
field (1995) and core convection in rotating A-type 
stars (Browning et al. 2004). 

The work by Freytag et al. (1996) motivated 
Herwig et al. (1997) to include the exponential 
overshooting derived from hydrodynamical simu- 
lations to all convective boundaries. The e-folding 
distance for convective boundaries in the stellar 
interior was derived semi-empirically for main- 
sequence stars, and then applied to AGB stars 
(Herwig 2000). Similar approaches were consid- 
ered by Mazzitelli et al. (1999), Ventura et al. 
(2000), Cristallo et al. (2004), and Althaus et al. 
(2005). 
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However, convection in shells driven by nuclear 
energy release is different than envelope convec- 
tion. Typically, the Mach numbers are smaller 
throughout the convection zone, and the border- 
ing radiative layers are relatively more stable. Ac- 
cordingly, the semi-empirical determination of the 
overshooting value /ov 1 gives smaller values for 
the stellar interior (/ov < 0.016) compared to 
the shallow surface convection (0.25 < /ov < 1-0) 
according to Freytag et al. (1996). For deep en- 
velope (or core) convection / ov can be expected 
to be considerably smaller since the ratio of the 
Brunt- Vaisala, time-scales of the stable to the un- 
stable layers decrease with increasing depth, i.e. 
adiabaticity (Herwig et al. 1997). 

Relevant in this regard are the 2D simulations 
of interior O-burning shell convection in massive 
stars immediately before a supernova explosion 
of Bazan & Arnett (1998) and Asida & Arnett 
(2000). In these simulations hydrodynamic fluid 
motions are not confined to the convectively unsta- 
ble region. The authors describe the excitation of 
gravity waves in the adjacent stable layers (Hurl- 
burt et al. 1986), as well as an exponential decay 
of convective velocities at and beyond the convec- 
tive border. As a result highly reactive fuel in the 
form of C is mixed from above into the O-shcll. 
However, the authors caution that the quantita- 
tive results may suffer some uncertainty due to 
numerical effects like resolution. The main dif- 
ference between these simulations and our regime 
of the He-shell flash is the lower driving energy 
in our case, and we expect correspondingly lower 
Mach numbers, less mixing and less overshooting. 
Young et al. (2005) present initial results on up- 
dated O-shcll 2D and 3D convection simulations 
and apply their findings on hydrodynamics mix- 
ing to ID stellar evolution models of supernova 
progenitors (Young & Arnett 2005). The inter- 
nal structure in these models is markedly different 
from standard models, with important implication 
for the supernova explosion properties. 

1 Here, the overshoot formulation is Oov = 
Do cxp ^ f~ 2 j] j i where Do is the mixing- length the- 
ory mixing coefficient at the base of the convection zone, z 
is the geometric distance to the convective boundary, H p is 
the pressure scale height at the convective boundary, and 
/ov is the overshooting parameter (Herwig et al. 1997). 



1.3. About this paper 

The main purpose of our study is to obtain a 
characterization of the hydrodynamic properties 
of shell convection in low-mass stars which is de- 
termined by weaker driving due to smaller nuclear 
energy generation compared to massive stars. We 
investigate the typical morphology, the dominat- 
ing structures and the time-scales of shell-flash 
convection, and how these properties depend on 
resolution and driving energy. We identify areas 
of parameter space that require further study in 
two- and three-dimensions. In our present work 
we do not intend to quantify mixing across the 
convective boundary, but we want to gain insight 
for future simulations that shed more light on this 
important question. 

In the following section we describe the ID stel- 
lar evolution code and the hydrodynamics codes 
used in this work. Sect. 3 describes the stellar evo- 
lution sequence and the specific model that serves 
as a template for the hydrodynamic simulations. 
Then we describe the results in Sect. 4; specifically 
we describe the general properties of the convec- 
tion, the onset of convection, the evolved convec- 
tion and the convective and oscillation modes ob- 
served in our simulations. We close the paper with 
a discussion (Sect. 5) of the results, in particular 
the comparison of the ID MLT velocity profile and 
the time-averaged vertical velocity profile from our 
hydrodynamic simulations. 

2. Codes 

2.1. Stellar evolution code 

Our ID hydrostatic, Lagrangian stellar evolu- 
tion code EVOL with adaptive mesh refinement 
and time stepping is equipped with up-to-date in- 
put physics (Herwig 2004, and references therein). 
Mass loss is included according to the formula 
given by Blocker (1995) with a scaling factor 
t]bl =0.1. Nucleosynthesis is considered by solv- 
ing a sufficiently detailed nuclear network, with 
reaction rates from the NACRE compilation (An- 
gulo et al. 1999). We use OPAL opacities (Iglc- 
sias & Rogers 1996) and low-temperature opaci- 
ties from Alexander & Ferguson (1994). In the 
stellar evolution code convective energy transport 
is described by the mixing-length theory (Bohm- 
Vitense 1958) with a constant mixing length pa- 
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rameter of Omlt = 1.7 determined by calibrating 
a solar model with the parameters of the sun. For 
material mixing we solve a diffusion equation for 
each species with a convective diffusion coefficient 
as given by Langer et al. (1985). 

2.2. Hydrodynamics code 

For the hydrodynamic simulations we use the 
multi-dimensional radiation-hydrodynamics code 
RAGE (Radiation Adaptive Grid Eulerian) de- 
signed to model a variety of multi-material flows 
(Baltrusaitis et al. 1996). The conservation equa- 
tions for mass, momentum, and total energy are 
solved through a second-order, direct-Eulerian Go- 
dunov method on a finite volume mesh (Dendy & 
Gittings 2005). RAGE has been extensively tested 
on verification problems (Belle et al. 2005; Holmes 
et al. 1999; Hueckstaedt et al. 2005). The code 
has the capability for continuous adaptive-mesh 
refinement to increase spatial resolution in areas 
of interest with a minimal increase in computa- 
tional time. This capability is in particular useful 
for multi- fluid applications. We do not use this 
feature because all our runs have a single-fluid. 

3. Stellar evolution models of the He- 
shell flash and initial multi-dimensional 
setup 

We simulate the He-shell convection zone at a 
time point during a thermal pulse evolution imme- 
diately before the peak of the He-luminosity. At 
this time the He-burning luminosity that drives 
the convection is already large implying a rapid 
build-up of an entropy excess that drives convec- 
tion. At the same time the layer has not yet been 
expanded too much, and the number of pressure 
scale heights involved is not excessive. Neverthe- 
less, our simulations vertically span 11.4H p . Our 
inital conditions are a a set of polytropic profiles 
that closely resembles the structure from the full 
ID stellar evolution model. 

3.1. The stellar evolution model 

Our hydrodynamic simulations are based on a 
model from run ET2 of Herwig & Austin (2004) 
that represents typical properties of a thermal pule 
in a low-mass asymptotic giant branchn star. This 
sequence has metallicity Z = 0.01 and the main- 
sequence initial mass is 2 M©. The evolution track 



was followed from the pre-main sequence until all 
envelope mass was lost at the tip of the AGB. We 
assumed exponential, time- and depth-dependent 
overshooting at the bottom of the convective enve- 
lope (/ov = 0.016). At the bottom of the He-shell 
flash convection zone a very small overshooting 
(/ov = 0.001) was applied, mainly for numerical 
reasons. Such a small overshooting value has lit- 
tle noticeable effect on the evolution but improves 
the convergence of the code. The diffusion coef- 
ficient drops 5 orders of magnitude within only 
7 grid points corresponding to about 7km, or less 
than O.OlHp. Sequence ET2 has 16 thermal pulses, 
eight of which are followed by a third dredge-up. 

We choose a model from the second to last 
(15 th ) thermal pulse of the ET2 run. At this ther- 
mal pulse, the core mass is 0.6 M© and the en- 
tire stellar mass is 1.4 M©. For a stellar model as 
template for the initial stratification of the hydro- 
dynamic simulations, we favor a large luminosity 
which implies large convective velocities. How- 
ever, the large luminosity quickly leads to a sig- 
nificant expansion of the layers above the He-shell 
(Fig. 1), which dramatically increases the number 
of pressure scale heights that the simulation would 
have to cover (Fig. 2). In particular, we would like 
to include enough stable layer above and below the 
convection zone to get an impression of the hydro- 
dynamic fluid behavior at the convective bound- 
aries and the adjacent stable regions. 

We picked model 70238 of that rerun as a tem- 
plate for our hydrodynamic initial setup. This 
model is one time-step (0.07 yr) before the peak 
of the He-shell flash luminosity of 4.29 • 10 7 L© 
(model 70239). The thermodynamic and mixing 
properties of model 70238 are given in Table 1 for 
the top and bottom of the convection zone and the 
location of where we chose the top and bottom of 
the hydrodynamic simulation box to be. The sim- 
ulations covers a total 11 H p , roughly half of this 
inside the convection zone. 

The molecular weight profile and the convec- 
tive velocity w CO nv = 5/(|aMLT^p) are shown 
for two consecutive models (70237 and 70238) 
in Fig. 3 and 4. The time step between the 
two models is 0.04 yr corresponding to approx- 
imately 1400 convective turnovers. The radius 
integrated He-burning energy in the two subse- 
quent models 70237 and 70238 is 4.7 • 10 17 and 
1.1 • 10 18 cmerg/g/s. During the He-shell flash, the 
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nuclear energy is mainly generated by the triple-a 
reaction. 

3.2. The hydrodynamic simulation setup 
and run parameters 

We simulate a plane-parallel layer in Cartesian 
coordinates which represents the intershell of the 
AGB star at the time nearly at the peak He-shell 
flash luminosity. In our models (IcO for 2D and lcl 
for 3D), the grid is equidistant in all dimensions, 
with a cell size of 55km in our standard mesh g 
(600x200). The box covers 33000km horizontally 
and 11000km vertically in the 2D runs. 

The bottom of the simulation box is located 
at a stellar model radius of 7507.5km, and subse- 
quently all length scales are with reference to the 
box bottom. The bottom of the convection zone 
is at 1650km. We approximate the stellar evolu- 
tion model stratification in the simulation box by 
a combination of three polytropic stratifications 
(Hurlburt et al. 1986), with 7 in each layer ap- 
proximating the 7 in the corresponding layer in 
the stellar evolution model. For the region be- 
low the convection zone we assume a stable poly- 
tropic stratification with 7 = 1.2. For the convec- 
tive region we set 7 = 7 a d = 5/3, while the top 
polytropic stable stratification mimics the stellar 
model with 7 = 1.01. At the bottom of the con- 
vection zone we set the temperature to 2.48 • 10 8 K 
and the density to 1.174-10 4 g/cm 3 as in the stellar 
model(see Fig. 6). We fix the temperature below 
the convection zone and at the top of the lower 
stable layer at 1.17 • 10 8 K, corresponding to a fac- 
tor of 0.47 between the peak temperature at the 
bottom of the convection zone and the local tem- 
perature minimum just below the convection zone. 
The location of the top of the convection zone 
is set to 7732.4km from the bottom of the box, 
corresponding exactly to the value in the stellar 
model. The top of the box at 11000km from the 
box bottom is chosen to be just below the H-rich 
envelope. Fig. 6 demonstrates the agreement be- 
tween our piecewise polytropic stratification and 
the stellar evolution template model. 

In the RAGE code the molecular weight of a 
material is specified in terms of the specific heat, 
which for an ideal gas is cy = 1 with the 
universal gas constant 1Z = 8.31 • 10 7 crg/K/g. 
fi— gradients are present at the bottom and top of 



the convection zone (Fig. 3). Although they add 
stability and reduce mixing across the boundaries, 
we ignore ^—gradients in this set of calculations. 
We plan to address the effect of /i-gradicnts in a 
future study. Here, we chose as a single mean 
molecular weight fi = 1.4 which closely matches 
the value in the convection zone from the stellar 
model. Finally, we assume a constant gravity of 
log g = 7.7[cgs] vertically directed downward as in 
the convection zone of the stellar model. In the 
stellar model gravity varies from log g = 8.05 at 
the bottom of the simulated layer to logg = 7.35 
at the top. 

In the stellar evolution model, the nuclear 
energy that drives the convection is calculated 
by solving a nuclear network based on realistic 
temperature-dependent nuclear reaction rate in- 
put. The same approach used in the hydrody- 
namic calculation may introduce additional res- 
olution dependencies of localized flow morpholo- 
gies, in particular at the bottom of the convection 
zone. We assume here a time-, temperature- and 
resolution-independent volume heating in a small 
layer at the bottom of the convection zone that 
releases the same integrated amount of energy as 
the stellar evolution model. Specifically, we release 
approximately e = 2 • 10 10 erg/g/s over 550km at 
the bottom of the adiabatic stratification. A con- 
stant volume heating is realized by multiplying eo 
with the density at the bottom of the convection 
zone for the input value (e* = eopo) and then 
locally in each cell adding e = e*/p to the cell 
energy. 

We introduce random density and temperature 
perturbations in pressure equilibrium on the grid- 
level of logAT/T = —4 and centered at the un- 
perturbed values. Due to the discretization er- 
ror, the initial numerical stratification slightly vi- 
olates hydrostatic equilibrium. This would cause 
a spectrum of plane-parallel p-mode oscillations 
with significant velocity amplitudes. During an 
initial damping phase (700s) of the simulations we 
damp those artificial p-modes to amplitudes at the 
cm/s- level. The damping is achieved by multiply- 
ing a factor to the momentum in each cell. This 
damping factor decreases from an initial value to 
zero over the damping phase. A corresponding- 
heating factor increases from to 1 over the damp- 
ing phase and smoothly turns heating for convec- 
tion driving on. The length of the damping phase 
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was chosen so that an overall velocity minimum 
could be achieved before the onset of convective 
motions. 

We formulated a grid of runs in 2D that cover 
a range in resolution and heating rate (Fig. 5). 
These runs are labeled IcOAB, where lc = lumi- 
nosity driven convection, is the run series, in 
this case a series of 2D runs, A denotes the heat- 
ing rate and B the grid size. We have done one 
3D run, denoted by lcl. We also discuss two runs 
(lc2) in which we have used three formally differ- 
ent species in the three layers of our setup, but the 
three materials have the same material properties. 
This allows some preliminary insight into mixing. 

3.3. Time-scales and summary 

The convective turn-over time scale is of the 
order of 600s and convective velocities are of the 
order a few km/s. We use an explicit, compress- 
ible hydrodynamics code, in which the time step 
is limited by the Courant-Fricdrichs-Levy condi- 
tion imposed by the sound speed. The maximum 
sound speed in our setup is > lOOOkm/s. Accord- 
ingly, our convective flows have a low Mach num- 
ber of the order 10 _3 Ma. With a grid size of 55km 
this implies a time step of At < 5 • 10~ 2 s. In or- 
der to achieve some convective steady-state and 
evolve this configuration for some time to have 
good statistics, our multi-dimensional simulations 
run for approximately 15 to 30 convective turn- 
over times, which corresponds to more than 10 5 
time steps. Thus, the multi-dimensional simula- 
tion time is much less than 0.1% of the stellar evo- 
lution calculation time-step that reflects the ther- 
mal time scale. Our calculation represent a snap- 
shot in time compared to the long-term evolution 
of the entire Hc-shcll flash which takes 200—300 yr. 
In our simulations we assume a constant heating 
rate in time. For the highly adiabatic He-shell 
flash convection the ratio between convective and 
non-convective heat transport (the Peclet number) 
is very large. For that reason we ignore thermal 
transport in this set of simulations. 

By running a series of simulations with heating 
rates spanning a factor 1000 we study how hydro- 
dynamic properties depend on the heating rate. 
We also want to derive how our results depend on 
resolution. Fig. 5 gives an overview of our 2D runs 
for our standard single-fluid setup IcO. 



4. Results 

Our simulations start with a Raleigh- Taylor like 
growth sequence of the initial perturbations driven 
by the heating at the bottom of the unstable layer 
(Sect. 4.1). Eventually, fully developed convection 
emerges (Sect. 4.3). It features a cellular flow con- 
figuration best seen in the pseudo-streamline plot 
in Fig. 7. The flow field is dominated by a few large 
cells (about 4 in that example) each consisting of 
a clockwise and a counter-clockwise vortex. Verti- 
cally the cells span the entire convectively unsta- 
ble region. However, the vortex centers are located 
closer to the bottom, and leave room for some sub- 
structure in the top layers of the convectively un- 
stable zone. Additional smaller rolls also appear 
close to the bottom boundary. Gravity waves are 
excited in the stable layers. An oscillaton mode 
analysis of the their appearance at the transition 
from unstable to stable stratification can be used 
to better understand the interaction of stable and 
unstable layers (Sect. 4.5). Fig. 8 shows the corre- 
sponding entropy fluctuations for the same snap- 
shot. Hot material rises from the heating layer 
at the bottom of the convection zone. The bulk 
of the high-entropy areas starts to turn side-ways 
well below the top of the unstable layer. However 
vertical high-velocity streams persist through the 
unstable layer, and the velocity field can be per- 
pendicular to the streaks of high-entropy material. 

4.1. Onset of convection 

The drag force applied in the first 700 s (p- 
modc filter, Sect. 3.2) of a simulation eliminates 
efficiently all modes that appear due to the imper- 
fect hydrostatic equilibrium of the initial stratifi- 
cation. On the other hand, it also delays the onset 
of convection somewhat. 

After some time, the external heating builds 
up a small extra entropy bump at the bottom 
of the entropy plateau. The top of this bump 
is the place where the entropy has the steepest 
negative slope, and where the convective instabil- 
ity causes small-scale modes to grow first (see the 
top panel in Fig. 9). The modes expand verti- 
cally and the velocity amplitude increases. Later, 
when they leave the linear regime, they form a 
patch of tiny mushroom-like instabilities that are 
similar in their morphology to the Raleigh- Taylor 
instability. These grow and merge further forming 
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complex large-scale structures (Fig. 9). 

The stable regions are far from being quiet but 
instead filled with waves. Remarkably, these waves 
occur in the upper stable zone well before the 
convective plumes reach this layer. The motions 
are not excited by plumes "hitting" the bound- 
ary layer but by the pressure excess in front of 
the plumes. This can be clearly seen in Fig. 10 
which shows the pressure fluctuations and pseudo- 
streamlines of the same time and model as the 
middle panel of Fig. 9. The rising plumes have not 
yet reached the top of the unstable layers; however 
the entropy fluctuations at and above the top of 
the unstable layer show that oscillation modes are 
already excited by the fluctuating over-pressure 
regions seen, for example, at the horizontal po- 
sition x—25 Mm in Fig. 10. We demonstrate in 
Sect. 4.5 that these oscillation modes are in fact 
gravity waves. 

Figure 12 demonstrates that the amount of de- 
tail visible in the plumes strongly depends on the 
resolution while the number of plumes (or the 
size of the largest structures at a given time-step) 
is only slightly affected. The wavelength of the 
modes that initially start to grow and also the 
size of the first generation of mushrooms decreases 
with grid size. However, these small-scale struc- 
tures quickly merge to form much larger cells. 

We do not present models with a coarser grid 
than 210x70 because test runs have shown that 
numerical problems arise at the upper boundary 
of the box (an artificial entropy inversion and tiny 
"surface cells"). That is not particular worrisom; 
a hydrodynamics code needs a minimum number 
of grid cells to resolve a pressure scale height prop- 
erly, and we have about 11.4H p inside our box. 

In the low-resolution model (top panel in 
Fig. 12) the entropy above the "mushrooms" is 
not constant as it should be for perfectly adia- 
batic flows. This artifact essentially vanishes with 
increasing resolution and docs not seem to have 
any adverse effect. 

To demonstrate the influence of the heating 
rate upon the convective patterns during the onset 
phase, a scaling of the amplitude of the entropy in- 
homogeneities (the fluctuations increase with the 
heating) and adjustment of the selected time-steps 
are necessary (a larger heating rate causes larger 
convective velocities and shorter growth times, 



see Fig. 13). Taking the scaling into account, 
the resulting morphology of the convective flows 
are remarkably similar, suggesting the possibility 
of avoiding excessively CPU-intensive simulations 
with low heating rates and correspondingly small 
convective velocities. 

Careful inspection of the lower panels of Fig. 13 
reveals spurious entropy minima as dark spots at 
the top of the entropy plateau caused by a small 
negative artificial entropy spike at that position. 
While the employed initial p-mode filter facili- 
tates the relaxation to a complete numerical hy- 
drostatic equilibrium, there is also a thermal re- 
laxation going on with a much longer time-scale 
(hundreds of seconds instead of seconds). This 
process leads to a smearing of the initially sharp 
jump in density and internal energy at the bottom 
of the entropy plateau and in the discontinuity in 
the derivatives at the top of the plateau. This 
redistribution of density and energy does not pre- 
serve the entropy profile but causes an artificial 
(small) entropy maximum and minimum to occur 
at the bottom and the top of the entropy plateau, 
respectively. The amplitude and the width of the 
entropy spikes decrease with increasing resolution. 
However, at low resolutions and low (i.e. realistic) 
heating rates, they cause an initial driving of the 
convective flow comparable to the one due to the 
external heating. 

Fortunately, this additional driving is a tran- 
sient phenomenon because there is only a lim- 
ited amount of energy available that can be "ex- 
tracted" by a conservative code out of the discon- 
tinuities. It only increases the time it takes for 
the convection to settle to a statistically station- 
ary state by approximately 3000 sec (for the run 
IcOgg) depending on resolution. 

4.2. Onset of convection: 2D versus 3D 

In addition to the two-dimensional models, we 
have some information from a single 3D run (lcld, 
1502x100 grid points). Figs. 14 and 15 display 
entropy inhomogeneities during the onset of con- 
vection for the 2D and the 3D run, respectively. 
Heating rate, vertical extent, and grid cell size are 
the same in both models. 

The 3D model initially has smaller velocities 
and appears to be slightly behind the 2D model. 
However, as soon as convection sets in, the ve- 
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locities in the corresponding snapshots are very 
similar. The structures in the last panels seems 
to be slightly smaller in 3D than in 2D. However, 
distorted "mushrooms" appear in both sequences. 
Also, gravity waves are excited in both models. 

We do not sec any hint in the 3D models that 
would indicate that the results we derive from the 
grid of 2D models would suffer from the restriction 
to two dimensions. To get accurate quantitative 
results 3D models are to be preferred. However, 
3D models with high resolution and realistically 
low heating rates are computationally extremely 
demanding. We will try to approach this regime 
in future work. 

4.3. Evolved convection 

We demonstrate the properties of fully devel- 
oped convection by means of sequences of snap- 
shots and averaged ID profiles. 

Figure 16 shows an example of the evolution 
of the convective patterns in a later stage of the 
1200x400 run with realistic heating rate. During 
this time, the flow is dominated by 4 global cells 
(the aspect ratio of the unstable zone is about 6). 
The centers of the vortices are located in the lower 
half of this zone, characteristic for convection in a 
strongly stratified environment, as found by Hurl- 
burt et al. (1984) in their stationary convection. 
However, in our case the upflows, starting in the 
thin heating layer at the bottom, are dynamic and 
change their shapes all the time. Smaller rolls can 
occur, especially at both sides of the foot-points of 
the updrafts and also near the top of the entropy 
plateau. 

The stability of the general pattern is supported 
by the fact that the down-flows tend to inhibit 
the development of new updrafts that can develop 
due to the convective instability in the heated bot- 
tom layers. Instead, the heated material is pushed 
sideways into the existing upflow channels. Dur- 
ing this process, new upflow "fingers" can develop. 
Most often, they are washed into the larger up- 
drafts. However occasionally, they have a chance 
to grow and to change the global pattern. 

The situation is pretty much the inverse of 
the well-known scenario of stellar surface convec- 
tion, where radiative cooling creates a thin super- 
adiabatic layer that leads to the formation of cool 
downdrafts that drive the convection. Whereas in 



surface convection the stratification tends to focus 
the downdrafts (Nordlund et al. 1997), the hot up- 
flows expand. This expansing effect can widen the 
thin columns of hot upflowing material in the case 
of He-shell flash convection. Nevertheless, they 
can rise over several pressure scale heights before 
they break up and can traverse the entire convec- 
tion zone. 

The entropy jump at the bottom of the plateau 
and the sudden increase at the top (see Fig. 6) con- 
fine the convective flow essentially to the region in 
between - with very narrow transition layers. In 
these boundary layers the largest differences be- 
tween the runs with different resolutions and heat- 
ing rates occur. To compare them quantitatively, 
we compute some averaged quantities. 

During the simulations a snapshot is stored ev- 
ery 0.5 s. This results for each quantity q (the ve- 
locity components, internal energy, p, T, P) in a 
data-set q(x 1 y, t) depending on horizontal coordi- 
nate x, vertical coordinate y, and time t. This out- 
put high rate clearly over-samples the slow convec- 
tive motions especially in the runs with low heat- 
ing rates. However, it is justified for the fastest 
gravity waves. And it is just barely sufficient to 
sample the pressure waves (see Sect. 4.5). 

For each snapshot we compute horizontal aver- 
ages, 

9av g (y,£) = {q(x,y,t)) x (l) 

and horizontal root-mean-square (rms) fluctua- 
tions 

/ 2 \ 1/2 

qAi ms (y,t) = [{{q(x,y,t) - q avg (y,t)) } x j 

(2) 

We check the time evolution of these quantities to 
ensure that the simulations have reached a statis- 
tical steady state. Due to heating at the bottom of 
the convective layer the entropy of the convection 
zone is slowly increasing and the extent of the con- 
vection zone is growing at a small rate. The time 
evolution of the maximum vertical rms- velocity is 
shown for runs with standard and enhanced heat- 
ing and different resolutions in Fig. 11. For the 
runs with the standard heating rate g a steady- 
state according to this disgnostic is reached af- 
ter 6000s, following the onset phase of convection 
which is characterized by larger velocities. At a 
higher heating rate the difference between the ini- 
tial velocity peak during the onset phase and the 
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following steady-state velocity is smaller. In this 
case the steady state is already reached after ap- 
proximately 3000s. 

To compare models with different resolutions 
or heating rates we average these quantities over 
time to get time- independent vertical profiles, as 

q avg (y) = (<?av g (y,i))t ( 3 ) 

or 

1 /2 

<7Arms(?/) = ((<7Arms(y,i) 2 )t) • (4) 

Examples for the rms-fluctuations of vertical and 
horizontal velocity are shown in Fig. 17, first and 
second row, respectively. Plotting the entropy s 
on the same scale as in Fig. 6 would show no dif- 
ference to the initial stratification. Therefore, only 
the entropy difference to the minimum plateau 
value is plotted in the panels in the third row, on a 
highly magnified scale. To allow the comparison of 
models with very different heating rates (and en- 
tropy fluctuations) in the last row the same data is 
rescaled to give the entropy bump at the bottom 
of the plateau a value of unity. 

The velocities are clearly not confined to the re- 
gion that would have been predicted by MLT (be- 
tween 1.7 Mm and 7.7 Mm, compare the top pan- 
els in Fig. 17 with Fig. 4). Nevertheless, before 
reaching the boundaries of the entropy plateau 
the vertical velocities drop exponentially by orders 
of magnitude. Close to the boundaries there are 
strong rapidly declining horizontal velocities lead- 
ing to significant shear flows (second row of panels 
in Fig. 17). 

Although the classical boundaries of a con- 
vection zone do not actually confine the velocity 
field completely, they still indicate where the flow 
changes its character: large overturning convec- 
tive motions, inside and small wave-like motions 
outside (cf. Fig. 17 and e.g. Figs. 7 and 8). The 
low-heating models in Fig. 17 (top right panel) in- 
dicate the presence of a transition region where 
the velocities decline rapidly inside the top stable 
regions before they rise due to the waves. 

The dependence of the velocity on the heat- 
ing rate is rather simple: it grows monotonically 
with heating rate. The waves in the stable layers 
have larger amplitude than the velocities of the 
overturning convective motions (right column in 
Fig. 17). It is important to emphasize that the 
large vertical velocities shown in the stable lay- 



ers of the top panels of Fig. 17 do not necessarily 
indicate effective mixing (see Sect. 4.6 for details). 

The dependence on the resolution is somewhat 
similar. For the convective motions, there are 
hints towards convergence of the vertical and hor- 
izontal rms-velocities at the highest resolutions 
models (1200x400). However, the wave amplitude 
in the runs with large heating rate shows less con- 
vergence. The (logarithmic) step in the amplitude 
of the horizontal velocities between the 1200x400 
and the 600x200 run is smaller than between the 
600x200 and the 300x100 run. Models with even 
higher resolution are desirable to check the con- 
vergence. The appearence of similiar quantitative 
properties in He-shell flash convection appears to 
be robust; these properties are present at all reso- 
lutions and heating rates, in 2D and in 3D. 

A higher resolution does not affect the number 
of large cells noticeably. However, it leads to finer 
entropy filaments and complex sub-structures of 
the largest plumes (cf. Figs. 18 and 19). Less 
numerical dissipation on a given scale leads to 
higher velocities, stronger velocity gradients (caus- 
ing more dissipation) , and some additional smaller 
rolls. However, even the high-resolution model 
in Fig. 7 shows the dominance of the large-scale 
structures. The trend to higher convective veloc- 
ities is limited by the amount of energy that has 
to be transported. 

The situation is different for the velocity am- 
plitudes of the g-modes in the stable layers. A 
wave with a certain wavelength directly feels the 
lower damping rate at a higher numerical resolu- 
tion. Since the wave excitation rate is larger due to 
the larger convective velocities, a mode reaches a 
significantly larger amplitude. This explains why 
the g-mode velocity amplitudes increase with res- 
olution. 

4.4. Detailed shape of the entropy plateau 

The initially flat entropy profile on the plateau 
changes during a simulation. However, it does not 
show a simple monotonic decline as in a corre- 
sponding MLT model. Instead, it has a complex 
height and parameter dependence (see the panels 
in the lower half of Fig. 17). 

All models show a sharp entropy peak in the 
heating zone at the bottom of the entropy plateau. 
The entropy declines with height, indicating a con- 
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vectively unstable region consistent with MLT. 
However, the entropy minimum is not located at 
the top of the plateau but somewhere in the mid- 
dle, typically in the lower half. Above the mini- 
mum, the entropy rises and has, in most models, 
a smooth monotonic transition to the stable lay- 
ers above the convection zone. An exception is 
the model IcOgh (with realistic heating rate and 
highest resolution, cf. the two bottom panels in 
the center column in Fig. 17): It shows a second 
minimum at the top of the plateau. 

Note, that even if the entropy gradient is pos- 
itive in some parts of the convection zone and 
therefore the stratification is formally stable there, 
the matter is heated in the shallow zone at the 
bottom and the plumes maintain their entropy ex- 
cess and remain buoyant over the entire convection 
zone (indicated by the bright color of the plumes 
in e.g. Fig. 8). 

We interpret the entropy rise at the top of the 
plateau of most models as due to mixing of high- 
entropy material from the stable layers above into 
the convection zone. This mixing depends on the 
velocity amplitude, the numerical resolution, and 
time. 

A significant part of the flow close to the upper 
boundary of the entropy plateau is due to wave 
motions. These are almost reversible and have a 
low mixing efficiency and cause a diffusion of en- 
tropy due to numerical viscosity. In addition ro 
numerical mixing there may be a finite amount of 
physical mixing. 

While higher resolution means less dissipation 
it also means larger amplitude of the motions that 
cause the mixing. A large heating rate means 
strong mixing but also a quick rise of the plateau 
so that the relative effect remains small. The 
model IcOcg (which has the largest heating rate) 
has the largest absolute amount of mixing (the 
highest entropy rise at upper part of the plateau, 
see third panel in right column in Fig. 17). How- 
ever, the relative effect is smaller than for the mod- 
els with lower heating rate (fourth panel in right 
column in Fig. 17). 

The resolution effect onto the models with large 
heating rate (left column in Fig. 17) is only mod- 
erate. The relatively quick rise of the entropy 
plateau value seems to limit the region with the 
entropy contamination from above. The dissipa- 



tion of entropy is similar for all four tested resolu- 
tions. 

In the models with realistic heating rate 
(center column in Fig. 17) the entropy profile 
changes qualitatively and quantitatively. The low- 
resolution models show a strong rise of entropy, 
whereas the highest resolution model has a well- 
defined minimum close to the top of the plateau. 
The effect of this minimum can be seen e.g. in 
Fig. 8 or Fig. 16 as tiny dark plumes emanating 
from the top unstable layers downward. An initial 
extra minimum in the entropy profile of model 
IcOgg (600x200 points) vanishes during the simu- 
lation. The minimum of model IcOgh (1200x400) 
remains during the available simulation time. The 
profile for very high resolutions and long simula- 
tion time is therefore not predictable, yet. 

The change in the entropy profile due to this 
slow mixing has an impact on the flow for the runs 
with low resolution and low heating rate. While 
the flow patterns during the onset and early stages 
of convection are very similar when a proper scal- 
ing of the amplitudes of the fluctuations is taken 
into account (see Fig. 13) the patterns differ after a 
sufficiently long simulation time (see Figs. 20 and 
21). The contamination of the convection zone 
with high-entropy material is so strong that the 
stratification is partly stabilized. The convective 
velocities in the top part of the entropy plateau 
decrease, and the typical cell sizes become notice- 
ably smaller. 

4.5. Modes 

To further investigate the waves encountered in 
the simulations, we construct k-u> diagrams from 
Fourier analyses. An example is shown in Fig. 22 
for an arbitrary height point inside the upper sta- 
ble region for the run IcOgg. 

In order to produce such a diagram, we col- 
lect all points for a single quantity (the vertical 
velocity) for a chosen height and all time steps 
of an interval into a 2D array. For the example 
from the IcOgg run, the array has a space axis 
with 600 points spanning 33 Mm, and a time axis 
with 12000 points spanning 6000 s. Next, a pos- 
sible trend in time is removed by subtracting a 
parabola fitted to the time-dependent spatial av- 
erage of the data. This is actually not necessary 
for the vertical velocities. But it is essential for 
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data like the entropy that has a large average value 
and a trend in time. To avoid signals due to the 
misfit of the beginning and end of the sequence, 
the data is faded in and out by multiplying with a 
cosine-bell function. A 2D Fourier analysis is ap- 
plied and the power-spectrum is computed. The 
latter is smoothed in frequency space by apply- 
ing a (1/4, 1/2, 1/4) filter. To emphasize small 
signals, the power -0 - 1 is plotted in the examples. 

A simple mode with a certain wavelength and 
period would give a single (dark) spot in a k-uo di- 
agram. Imperfections (e.g. due to interaction with 
convection or damping) would cause a smearing in 
frequency. Plane-parallel modes are found on the 
v axis. Static signals (e.g. processes on time-scales 
exceeding the interval) are located on the 1/ A axis. 

While a horizontally moving sine wave would 
give a single spot, a more complex function that 
travels horizontally with a certain speed would 
produce an entire ray with a slope v\. Mode fam- 
ilies give rise to ridges. Their detailed shape de- 
pends on the background stratification. Convec- 
tion does not consist of discrete modes. Instead, it 
has only typical time and length scales and gives 
a smeared blob. 

On the left side of Fig. 22, pressure waves (p- 
modes) are visible, reaching up to the top of the 
diagram. In fact, a reflection of the pattern oc- 
curs at the top of the figure. The sampling rate is 
not quite high enough to fully resolve all p-modes, 
which leads to this aliasing phenomenon. 

The bottom part of the diagram (below the 
Brunt- Vaisala-frequency of 0.1s -1 ) is dominated 
by gravity waves (g- modes). At longer horizon- 
tal wavelengths (small 1/A) individual ridges are 
visible. The top ridge (highest frequency) corre- 
sponds to the mode with largest vertical wave- 
length. Modes with smaller wavelength have a 
lower frequency (contrary to the behavior of p- 
modes). The g-modes with very small horizontal 
or vertical wavelength are not resolved individu- 
ally and only give a dark "band" . 

The inconspicuous convective signature is the 
inclined flattened streak in the very lower left cor- 
ner at A>lMm, i>100s. To exhibit the location 
of the convective signal in the k-u diagrams more 
clearly Fig. 23 shows only the lower left corner of 
the entire diagram, with a non-linear scale that 
magnifies the region corresponding to long time- 



scales and large wavelength even further. The fig- 
ure shows a sequence of k-u; diagrams for various 
heights inside the stable regions and the convec- 
tion zone. 

All plots show p-modes with varying relative 
power. Inside the convection zone (y=4.70Mm) 
they are barely visible above the huge background 
of the convective motions. However, in the stable 
layers close to the top of the computational do- 
main they have almost the same total power as 
the gravity waves. The top panels in Fig. 23 show 
a single ridge between p-modes and g-modes: the 
fundamental mode (f-mode). 

The gravity waves in the lower cavity (four 
top panels) have a higher Brunt- Vaisala-frequency 
than the modes in the upper cavity. Both wave 
families can be traced into the convection zone 
(y=1.79Mm and y=7.45Mm). 

The extended convective blob is visible in all 
panels, reaching from the convection zone, where 
it dominates the power, far into the stable regions. 
A k-u diagram even closer to the top than the last 
panel in Fig. 23 does not show a clear signature of 
the convective velocities anymore (the location of 
the last panel was chosen accordingly). 

While the convective flow overshoots into the 
stable regions, the stable regions "shoot back", 
they tunnel somewhat into the region with the 
essentially flat entropy profile. This overlap of 
flows leads to the mixing of high-entropy material 
into the convection zone, as described in Sect. 4.3. 
The region where convective flows and gravity 
waves can interact is vertically extended and not 
confined just to the top layer of the convection 
zone. During the onset phase of convection, grav- 
ity waves are excited in the upper stable region 
well before the convective upflows reach the up- 
per boundary of the entropy plateau. Convective 
elements (plumes) do not just "smash" against 
the boundary. Instead, the braking of the flow 
(and the excitation of the waves) occurs over some 
extended region. Rising co-moving fluid pack- 
ets (adiabatically) compress the material between 
themselves and the top boundary. The pressure 
rises, brakes the vertical movement of the convec- 
tive element, and accelerates material horizontally. 
As a result, the vertical motion of the convective 
element is turned into a horizontal flow and the 
over-pressure at the top lifts the boundary of the 
entropy plateau slightly. The material in the sta- 
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ble region is pushed away both horizontally and 
vertically. The reaction of the stable layers has an 
oscillatory component - the g-modes. However, as 
long as the local over-pressure from the convection 
zone below exists, some distortion remains in the 
stable layers. The distortion closely follows the 
time-evolution of the convective flow and has the 
same signature in a k-u diagram. 

While the picture as described above remains 
qualitatively the same for all resolutions and heat- 
ing rates in our grid, there arc some differences and 
trends visible. The position of a mode peak in a k- 
uj diagram remains essentially the same for all runs 
since it only depends on the background strati- 
fication which only changes slightly towards the 
end of the IcOcg run (largest heating rate). How- 
ever, more heating means higher convective veloc- 
ities and smaller time scales. Therefore, the rela- 
tion between time-scales of convection and modes 
changes. The convective blob moves closer to the 
ridges of the g-modes, and a separation becomes 
more difficult. 

The mode peaks appear sharper at lower heat- 
ing rates, where less interaction with convection 
and less change in the background stratification 
lead to longer mode lifetimes and "cleaner" peaks. 
Convective scales are slightly smaller at lower 
heating rates, indicated by a shifted maximum and 
stronger tail towards higher frequencies. 

An increase in resolution does not change the 
frequencies of already visible modes noticeably. 
However, more and more ridges appear and this 
indicates that the wide and flat g-mode nodes are 
not well enough resolved. 

4.6. Mixing 

Our mixing analysis is preliminary. Mixing in- 
side the convective region is directly correlated 
with the averaged vertical velocities (top panels 
Fig. 17). An important feature of these velocity 
profiles is the substantial decline of the velocities 
within the convective region, starting at a position 
several hundred km above the bottom and below 
the top of the convection layer. These plots in- 
dicate that the decay of the vertical velocity field 
is exponential. Outside the convective region, the 
velocity amplitudes of the g-modes dominate the 
averaged vertical velocity profile. Due to the wave 
nature of g-mode oscillations, these rather large 



velocities correspond to mixing that is orders of 
magnitude less efficient than that inside the con- 
vection zone, and possibly negligible far away from 
the convection zone. The fact that the g-mode ve- 
locity amplitudes superimpose onto the decay of 
the convective velocity field complicates deriving 
the effective mixing in and into the stable layers. 
A preliminary analysis of the oscillation modes 
(Sect. 4.5) suggests that the convective velocity 
field and therefor the corresponding mixing con- 
tinues the exponential decay into the stable layer. 
However, future analysis has to reveal the extent 
of the overshooting plumes interact with the hori- 
zontal motions of the g-modes, potentially leading 
to increased mixing in the stable layers beyond the 
convective boundary. 

We have performed two test runs with enhanced 
heating rate and low to moderate resolution (lc2df 
and lc2dg, Fig. 5) and with three individual fluids 
in the three different layers of the setup. The flu- 
ids are treated separately by the code but have the 
same properties, i.e. the same molecular weight. 
We use these fluids as passive tracer particles to 
establish an approximate upper limit on mixing 
(Hurlburt et al. 1994). We find that there is a 
finite mixing of material across both boundaries. 
In both runs at both boundaries there are typi- 
cally four horizontal layers involved in the main 
abundance transition. Mixing from the top stable 
layer into the convection zone seems to be more 
efficient for the higher resolution run. The top 
layer abundance reaches a mass fraction of approx- 
imately 10~ 5 inside the convection zone, whereas 
the lower resolved simulation (lc2df) reaches a few 
10~ 6 . Mixing of material into the convection zone 
from below is at a similar level but shows a signifi- 
cant upward trend at late times (after several con- 
vective turn-over times). Mixing of material from 
below the convection zone is roughly at the same 
level as mixing from above, but more efficient for 
the lower resolution case than for the higher reso- 
lution case. 

Mixing at the top unstable-stable interface can 
be seen in the model IcOdh with large heating 
rate and high resolution (cf. the bottom panel in 
Fig. 18). The shear flows at the top of the convec- 
tion zone "peel off" high-entropy material from 
the stable layers above. This material keeps its 
identity for some time, remains hotter than the 
surroundings, and is visible as bright spots in the 
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center of some vortices. The rather long survival 
time of these vortices may be a 2D artifact. 

The material in the lower stable zone is too 
dense to be mixed directly into the convection 
zone. However, the entropy (and density) jump 
is smoothed out over time. Then the flow at the 
bottom of some plumes becomes strong enough to 
lift up slightly cooler material from below that is 
finally mixed into the convection zone. It is visible 
as the dark areas directly beneath the plumes in 
Fig. 18. 

To some degree models with large heating rate 
and low resolution seem to have similar properties 
as models with smaller heating rate and higher 
resolution. We will therefore continue to perform 
simulations with too large heating rates in addi- 
tion to the runs with more realistic parameters. 

5. Discussion and conclusions 

We presented the first set of hydrodynamic sim- 
ulations specifically addressing the conditions in 
the convection zone close the peak luminosity of 
the He-shell flash. Compared to previous calcu- 
lations of other stellar convection zones our sim- 
ulations reveal some important differences. One 
of the most important is the stiffness of the con- 
vective boundaries of He-shell flash convection in 
contrast to the shallow surface convection of A- 
stars and white dwarfs as studied by Freytag et al. 
(1996). Other examples include solar-type convec- 
tion (Dintrans et al. 2005) or much of the parame- 
ter range in relative stability between the unstable 
and stable layers studied by Rogers & Glatzmaier 
(2005). Even in the simulations of core convection 
in A-type stars (Browning et al. 2004) convective 
plumes can cross the convective boundary signif- 
icantly. In shallow surface convection zones the 
convective downdrafts easily cross the convection 
boundary and penetrate deep into the stable lay- 
ers beneath. No significant g-mode excitation has 
been detected by Freytag et al. (1996), which is 
probably due to the absence of a stiff resonance 
floor with which convective plumes and the pres- 
sure field interact. However, it seems that because 
of this particular feature of the shallow surface 
convection zone Freytag et al. (1996) were able 
to study the decay of the convective velocity field 
in isolation. In our He-shell flash convection we 
do see exponential decay of the convective veloc- 



ity field as well. However it already starts inside 
the convection zone, and at or beyond the border 
quickly intermingles with the effects of the gravity 
waves. 

The comparison to the new O-shell convection 
simulations by Maekin & Arnett (2006, and 2006b, 
in prep) confirms the importance of g-modes for 
mixing and to correctly account for the physics 
at the convective boundary. Similar to our simu- 
lations their convective boundaries are very stiff, 
with overshooting of convective plumes into the 
stable layer reduced to a minimum. More impor- 
tant are horizontal wave motions which eventually 
induce turbulent mixing. A partial mixing zone at 
the convective boundaries is minimal. The O-shell 
convection simulations include ^-gradients, con- 
trary to our similations. [i gradients increase the 
impenetrable character of the convective bound- 
aries significantly, and we plan to include this ef- 
fect in our simulations in the future. 

We tentatively compare the mixing length con- 
vective velocity profile from our ID stellar evo- 
lution template model with the averaged vertical 
velocities in the hydrodynamic simulation with re- 
alistic heating rate in Fig. 24. Obviously, there are 
many differences in the assumptions of the two cal- 
culations, including geometry, micro- physics and 
of course the treatment of convection itself. Nev- 
ertheless, we find that the absolute value of the 
velocities inside the convection zone are similar. 
We do observe a different profile in the two rep- 
resentations. The hydrodynamic simulation show 
a more pronounced slope in the velocity profile 
within the convection zone than the MLT profile. 
The peak in the velocity profile from the hydro- 
dynamic simulation is close to the bottom of the 
convection zone. The most important difference is 
the significant decline of the velocity field already 
inside the convection zone. Near the lower bound- 
ary the vertical velocity is almost 3 orders of mag- 
nitude lower than the peak value. This may have 
implications for s-process nucleosynthesis as nuclei 
are on average exposed much longer to the hottest 
temperature at the very bottom of the convection 
zone. This should be considered in particular for 
cases that involve branchings sensitive to the con- 
vective time-scale of He-shell flash convection, like 
128 Xe (Reifarth et al. 2004). 

Consistent with previous work on stellar inte- 
rior convection we observe a rich spectrum of g- 
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modes excited in the stable layers. These g-modcs 
interact with the decaying convective flows at the 
boundary of the unstable region. Our simula- 
tions show that He-shell flash convection displays 
the main ingredients observed in previous stel- 
lar convection simulations, but with different rela- 
tive contributions of the various components. The 
actual overshooting of convective plumes across 
the convective boundaries is much smaller than in 
shallow surface convection and O-shell convection, 
while the excitation of internal gravity is probably 
stronger because the boundaries are stiffer. Al- 
though there is a significant overlap of the two 
processes we have shown that the analysis of the 
oscillation modes by means of the k-uu diagram al- 
lows us to disentangle the different contributions 
to the velocities. With this approach we observe 
that there is a finite mixing of convective motions 
accross the convective boundary. We postpone 
quantifying this mixing to future work. In ad- 
dition, we plan to quantify the contribution from 
wave mixing at He-shell flash convection bound- 
aries. 

We have identified a number of numerical prob- 
lems at low heating rates and low resolution that 
need to be dealt with in future 3D runs. These 
include the diffusion of entropy from above into 
the top zones of the unstable layer, impeding the 
development of the correct hydrodynamic flow- 
pattern. Our convergence and heating rate study 
shows that simulations with larger heating rate 
are in some ways equivalent to runs with higher 
resolution. Properties like the exponential decay 
rate of the convective velocity field seem to be only 
weakly dependent on the heating rate. However, 
the excitation of g-modes and their amplitudes de- 
pends sensitively on both the resolution and the 
heating rate. While we do see signs of convergence 
for the convective properties of our simulations, 
this is not the case for the g-modes. 

We consider 2D vs. 3D systematics, the effect 
of /i-gradients, a realistic treatment of nuclear en- 
ergy generation, and a more detailed study of the 
interaction of wave mixing with convective over- 
shooting as our most immediate tasks for future 
improvements. In addition we can make our sim- 
ulations more realistic by adopting an apppropri- 
ate stellar equations of state, by accounting for the 
spherical geometry of the He-shell, and by includ- 
ing a realistic gravitational potential. 
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Table 1 

Properties of He-shell region in stellar model 70238 





box bot. 


conv. bot. 


conv. top 


box top 


D/(cgs) 


0.000E+00 


1.634E+12 


1.051E+12 


0.000E+00 


R/R Q 


1.088E-02 


1.324E-02 


2.206E-02 


2.681E-02 


Hp/Ro 


1.296E-03 


2.439E-03 


1.258E-03 


2.337E-03 


M r /M Q 


0.5230 


0.5704 


0.5918 


0.5929 


p/(cgs) 


1.267E+05 


1.173E+04 


7.004E+02 


2.782E+01 


T/K 


1.544E+08 


2.479E+08 


4.631E+07 


4.872E+07 


lnP(cgs) 


48.68 


46.62 


42.16 


39.16 
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Fig. 1. — Time evolution of the radial location 
of the He-shell flash convection zone based on the 
ID stellar evolution model. Time is set to zero at 
the peak of the He-burning luminosity. The ordi- 
nate is zero at the stellar center. The grey shaded 
area represents the He-shell flash convection zone. 
Dots along the boundary indicate individual time 
steps in the ID evolution model sequence. The 
vertical solid line at t = — 0.07yr indicates the po- 
sition and extent of the hydrodynamic simulation 
box. Dashed and dotted lines correpsond to the 
Lagrangian coordinates given in the legend, and 
visualize the expansion of the He-shell as a result 
of the He-shell flash. 



Fig. 2. — Time evolution of the He-burning lu- 
minosity and the pressure at the boundary of the 
convection zone for the He-shell flash. The time 
axis, grey shades, and dots have the same meaning 
as in Fig. 1. 
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Fig. 3. — Mean molecular weight of the fully ion- 
ized material in the simulated shell of the two stel- 
lar evolution models before the peak He-burning 
flash luminosity. The radius is set to zero at the 
stellar radius of 7500km, coinciding with the bot- 
tom of the box for the hydrodynamic simulations. 
The two steps in the mean molecular weight are 
associated with the bottom of the He-shell flash 
convection zone at 1650km and the top at 7732km 
in model 70238. For model 70237 (dashed line) 
the bottom of the H-rich envelope is located at 
r = 10000km. 
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Fig. 4. — The velocity profiles from two subse- 
quent models of the ID stellar evolution model 
track using the mixing-length theory. Our hy- 
drodynamic simulations resemble the conditions 
of model 70238. 



Fig. 5. — The grid in resolution and heating rate 
of the standard 2D simulation set IcO. In the figure 
captions the combination of two letters identifies 
the simulation in this grid. The first letter gives 
the heating rate, the second the resolution. For 
example, run df has a heating rate of 30 times the 
standard case and a 300x100 grid. The numbers 
in the table give the length of the run in seconds. 
The standard heating rate g corresponds to what 
is described in Sect. 3.2. 
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Fig. 6. — The initial entropy, temperature, den- 
sity, and pressure stratifications of the hydrody- 
namic simulations (solid line) and the temper- 
ature, density and pressure from the ID model 
ET2/70238 (dashed line). 
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Fig. 7. — Fully developed convection in a high-resolution 2D run with standard heating rate. The flow 
field is represented by 25x75 pseudo-streamlines, integrated over 100 s for the constant velocity field of this 
single snapshot of evolved convection during the IcOgh run. The color indicates the corresponding pressure 
inhomogcneitics. The pressure itself does not deviate much from the initial values. Therefore, to make 
the fluctuations visible, the horizontal average of the pressure has been subtracted from the pressure value 
of every grid point. Bright means over-pressure (prominently where a "mushroom" approaches the upper 
boundary of the unstable region: compare with Fig. 8). Dark color indicates low-pressure (in the "eyes" 
of the vortices). The boundaries of the unstable layer (the entropy plateau in Fig. 6) at t/=1.7Mm and 
y=7.7Mm are clearly marked in the flow field and the pressure inhomogeneities. 



Fig. 8. — Entropy inhomogeneities for the same model of the IcOgh sequence as in Fig. 7. Again, the horizontal 
average of the entropy has been subtracted to render the small fluctuations visible. A bright color indicates 
material with higher entropy (and in fact higher temperature - due to the near pressure-equilibrium) . Dark 
means low entropy (or temperature). The boundaries of the entropy plateau at y=1.7Mm and t/=7.7Mm 
are again clearly visible in the change of the patterns. The subtraction of the horizontal mean causes bright 
features to be accompanied by dark horizontal stripes. These are pure artifacts of the visualization procedure 
and do not exist in the simulation data itself. 
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Fig. 9. — Entropy inhomogencities during the on- 
set of convection for the IcOgh run. 



Fig. 10. — Like Fig. 7 for an earlier time during the 
onset of convection. Note the build-up of higher 
pressure regions above rising plumes that can be 
associated with the excitation of gravity waves in 
the top stable layer. 



Fig. 11. — Time evolution of the maximum verti- 
cal, rms-averaged velocity for the standard heating 
rate g (top panel) and the enhanced heating rage 
d (bottom panel) for a range of resolutions. 



Fig. 12. — Resolution sequence of entropy inho- 
mogeneities during the onset of convection. 



Fig. 13. — Heating rate sequence (decreasing from 
top to bottom) of entropy inhomogeneities during 
the onset of convection. The time steps are scaled 
with q~ a - 21 . The color palettes are scaled with 
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Fig. 15. — 3D lcldf run sequence of the entropy 
inhomogeneities during the onset of convection. 
The same snapshot times and scaling is used as 
in Fig. 14. 



Fig. 14. — 2D IcOdf run sequence of entropy inho- 
mogeneities during the onset of convection. The 
same snapshot times and scaling is used as in 
Fig. 15. 



Fig. 16. — Time-sequence of the entropy fluctua- 
tions for evolved convection, IcOgh run. 
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Fig. 17. — Some time-averaged quantities for various model sequences. Each column of panels refers to the 
same sequence. All panels in one row show the same quantity. Left column: runs with different resolutions 
and large heating rate (IcOdX). The averages are performed over the interval 1500-4500s. Center column: 
runs different resolutions and realistic heating rate (IcOgX). The averages are performed over the (short) 
interval 3750-4300 s. Right column: runs with intermediate resolution (600x200) and different heating rates 
(IcOXg) . The averages are performed over the last 5000 s of each sequence. Top row: rms- fluctuations of the 
vertical velocity. Second row: rms-fluctuations of the horizontal velocity. Third row: mean entropy s shifted 
so that the minimum entropy on the plateau becomes zero. Bottom row: mean entropy s shifted as above 
and scaled to set the amplitude of the lower entropy bump to unity. 
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Fig. 18. — Resolution sequence and large heating 
rate, entropy inhomogeneities for evolved convec- 
tion. 



Fig. 19. — Resolution sequence and realistic heat- 
ing rate, entropy inhomogeneities for evolved con- 
vection. 



Fig. 20. — Heating rate sequence of the entropy in- 
homogeneities for evolved convection. The colors 
are scaled with q® 5 . 



Fig. 21. — Pressure inhomogeneities and pseudo- 
streamlines for evolved convection for runs with 
different heating rates. The colors are scaled 
with (jfj? 5 , the integration time for the pseudo- 
streamlines with 
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Fig. 22. — k-u> diagram (actually 1/X-i/ dia- 
gram) for an arbitrary height in the upper sta- 
ble region for run IcOgg. Both axes are linear 
and extend to the respective Nyquist frequencies. 
They are given by the sampling rate in time, 
^t,Nyqui S t=l/2Ai samp ii ng =ls _1 and the horizontal 
grid size, ^ i Nyquist=l/2Aa;g r id=9.09Mm _1 . The 
signal strength goes from bright (low amplitude) 
to dark (high amplitude). 
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Fig. 23. — k-uj diagrams for various heights of run IcOgg. The sample heights lie in the middle of the lower 
stable region, just below the entropy jump, inside and above the jump (top row from left to right); in the 
middle of the convection zone, near the top of convection zone, at the bottom of the upper stable zone, close 
to the middle of this zone (bottom row from left to right). 
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Fig. 24. — Comparison of vertical velocities de- 
rived from a hydrodynamics model (continuous 
lines: IcOgg) and mixing-length velocities from 
stellar evolution calculations (dashed lined: model 
70238). 
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